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Abstract 


Botero-Jaramillo, E., Romo, M. P., Méndez, B., €: Marengo, 
H. July-August, 2014). Use of the Distinct Element Method 
as a Tool to Detect Stability Problems in Deep Spillway 
Excavations. Water Technology and Sciences (in Spanish), 5(4), 
187-196. 


This paper describes a study to determine the stability of 
rock excavations in the La Yesca spillway channel using 
numerical analysis. The results of the analysis led to 
implementing important changes in the original conception 
of the spillway (reorientation, channel geometry and length) 
in order to ensure rock slope stability under static and seismic 
conditions. This paper describes the 3D numerical modeling 
of a 280 m high slope (corresponding to the original project) 
cut in a rock mass having complex geological formations. 
The modeling was carried out to evaluate both the static 
and dynamic stability of the projected rock slope. The 
construction site presented complex geological features such 
as faults and dikes as well as several geological formations 
resulting from varying qualities and fracture conditions. 
Available field data made it possible to incorporate the 
most important geological and topographical features. 
The spillway excavation was simulated using four stages. 
The effect of blast-induced damage and stress relaxation 
on the rock mass was also taken into account. After the 
simulation of the excavation process was completed, the 
seismic motion design for the La Yesca dam was input. These 
dynamic analysis simulated the three- dimensional behavior 
of the spillway which enabled identifying possible block 
instabilities, thereby determining whether to stabilize them 
or modify the original spillway layout. 


Keywords: Distinct element method, seismic analysis, slope 
stability, spillway problems. 


Introduction 


La Yesca is an important hydroelectric project 
(750 MW) located in north-western Mexico, 


Resumen 


Botero-Jaramillo, E., Romo, M. P., Méndez, B., € Marengo, H. 
(julio-agosto, 2014). El uso del método de los elementos distintos 
como una herramienta para la detección de problemas de estabilidad 
en las excavaciones del vertedero profundas. Tecnología y Ciencias 
del Agua, 5(4), 187-196. 


Este trabajo describe un estudio para determinar la estabilidad de las 
excavaciones en roca para el canal del vertedor de La Yesca utilizando 
análisis numéricos. Los resultados llevaron a la implementación 
de cambios importantes en la concepción original del vertedor 
(reorientación, geometría del canal y longitud), con el fin de 
garantizar la estabilidad de los taludes de roca en condiciones 
estática y sísmica. Se detalla el modelado numérico 3D de un 
talud de 280 m (correspondiente al proyecto original) construido 
en un macizo rocoso que tiene estructuras geológicas complejas. 
El modelado se llevó a cabo para evaluar la estabilidad estática 
y dinámica del talud proyectado en roca. El sitio de construcción 
presenta características geológicas complejas, como fallas y diques, 
así como varias formaciones geológicas resultantes de diferentes 
calidades de la roca y condiciones de fracturamiento. Los datos del 
campo disponibles hicieron posible la incorporación de la mayoría 
de las características geológicas y topográficas más importantes. La 
excavación del vertedor se simuló utilizando cuatro etapas. El efecto 
del daño inducido por las voladuras y la relajación de esfuerzos en 
la masa de roca también se tuvo en cuenta. Posterior a la simulación 
del proceso de excavación, se implementó el sismo de diseño de la 
presa La Yesca. Los análisis dinámicos simularon el comportamiento 
en tres dimensiones del vertedor, lo cual permitió la identificación de 
posibles bloques inestables, determinando de este modo si era posible 
estabilizarlos o modificar el diseño del vertedor original. 


Palabras clave: método de los elementos diferentes, análisis 
sísmico, estabilidad de taludes, problemas de vertedores. 


over the Santiago River, at the limit between the 
states of Jalisco and Nayarit. The construction 
of the spillway channel implied a cut as high as 
280 m height (see figure 1) in highly fractured, 
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geologically complex rock mass which was 
of great concern. The original design of the 
spillway consisted of a single channel 500 
m long and 96 m wide at its base. The main 
topographical and geological features of the 
site were captured in a three dimensional 
discrete element model that takes into account 
such features as the presence of rock fractures, 
faults, joints and geological contacts. The 
representative geological discontinuities were 
explicitly modeled and massive fracturing in 
the rock was accounted for by representing 
these materials as a rock mass, according to the 
Hoek-Brown criteria (1997). The poor hydraulic 
behavior of the spillway due to its partial 
obstruction due to falls of slope rocks possess 
a great risk to the dam safety. Accordingly, 
given the great importance spillways have 
on the safe behavior of dams, it was deemed 
of most importance to perform detailed 
numerical analysis to evaluate the slope safety 
of the high rock cuts. This article presents the 
analysis carried out for the original design of 
the spillway of La Yesca Dam and the decisions 
taken to decrease excavation-associated risks. 
The numerical model was first assessed 
by evaluating the stability of the rock mass 
before its excavation and then compared with 
the observed stable conditions of the rock 
mass over the time. Afterwards, the rock cut 


River flow 


Figure 1. Original project of the spillway channel. 


was simulated in four consecutive excavation 
stages. Stress states, displacement fields and 
spatial distribution of factors of safety were 
computed for each step. Finally, the dynamic 
stability analysis was carried out considering 
both the Operating Basis Earthquake (OBE) 
and Maximum Credible Earthquake (MCE), 
which are two synthetic motions proposed 
by the International Commission of Large 
Dams (ICOLD, 2010) representative of the 
site seismic environment (in this case of the 
Mexican subduction earthquake-generating 
mechanism) with 200 and 10 000 years of 
return period, respectively. On the bases of 
the analysis results, the risks associated to the 
spillway original lay out were assessed and 
decisions were taken. 


Site exploration 


The site geological materials characterization 
was carried out based on a broad geological 
and geophysical exploration (CFE, 2006). 
They are predominantly volcanic, but some 
alluvial and talus deposits were also found. 
The igneous rocks are heavily fractured and 
some dikes are present. The main lithological 
units include dacites, rhyolites, andesites, 
breccias and tuffs. Several faults are present in 
the region, being the Vertedor fault the more 
important. Figure 2 shows the fractures, faults 
and the spatial distribution of geomaterials at 
the main section of the spillway. The trace of 
the slope excavation is included as a reference, 
where the stepped line (blue line) depicts the 
initial excavation perimeter and the continuous 
line (red line) illustrates the average inclination 
of the model slope. The circled numbers allude 
to specific materials, which are presented in 
table 1. Given the length of the spillway, it is 
understood that not all of existing materials are 
visible in the section shown in this figure. 


Material properties 
Rock mass properties 


Triaxial tests were carried out on intact rock. 
The results were used along with the estimated 
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Figure 2. Geological structures at main section of the 


spillway. 


values of the Geological Strength Index (GSD), 
obtained from field observations, to determinate 
strength and deformability properties for the 
geomateriales using the Hoek and Brown 
failure criterion for jointed rock masses (Hoek 


was used for all materials and it was related 
with CoulomD's failure criterion to estimating 
an equivalent friction angle, Y”, and a cohesive 
strength, c', according to equations (1) to (3): 


k-1 
inp==— 1 
A a (1) 
0] 
C'= cm (2) 
2/k 
O, = 0,5" (3) 


Where k is the slope of the line relating the 
major and minor principal stresses 0/and 0”; O, 
and o, are the uniaxial compressive strengths 
of the rock mass and intact rock, respectively; s 
and q are constants which depend upon the rock 
mass Characteristics. The variation of q and c' 
were obtained as a function of o”, through a 


fitting process (Hoek éz Brown, 1997): 


€ Brown, 1997; Hoek, Carranza-Torres, « 
Corkum, 2002). The elasticity modulus E,, was 
estimated for two construction rock conditions: 
zero damage, D = 1 and D = 0.7, which is 
related to a moderate damage, induced by an 


A al 
60 m,(s+m, 9%,,) 


2(1+412+0)+6am,(s+m, Eo 


(4) 


sin $ = 


adequate and controlled excavation process 


(Brady $ Brown, 2004). The compression wave [d+20)5+(1-a)m,0,, ][s+m, o)” 


, 3n 3n 


velocity for the constitutive materials, V_, was e 0 


measured from geophysical explorations (see id o [sa m(s+m03,) | Jaja] 


table 1). The Hoek and Brown failure criterion (5) 


Table 1. Material properties. 


Material p (kg/m?) v mM, 

1. Dacite v, =3.8 km/s 2 500 0.35 6.25 Ss 
2. Dacite V, = 3.2 km/s 2.500 0.35 6.71 3 
3. Dacite Ye =2.7km/s 2 500 0.35 4.69 a 
4. Fractured breccia Ye. =2.5km/s 2500 0.35 3.34 Ss 
5. Uncompressed dacite 2500 0.35 2.37 Z 
6. Uncompressed breccia 2.500 0.35 1.78 + 
7. Dacite-rhyolite v, =1.7 km/s 2 500 0.35 1.64 3 
8. Dyke (andesite) 1 680 0.35 0.80 a 
9. Vertedor fault 2 500 0.35 0.93 S 
10. Vitreous tuff 2500 0.35 lo E 

= 

E 
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Where m, is the value of the constant m (see 
table 1) for the rock mass and o, = 0%, . /o,, is 
the upper limit of confining stress over which 
the relationship between the Hoek and Brown 
and the Mohr-Coulomb criteria are considered. 


3máx 


Its value was determined from the theory of 
elasticity, where the stress ratio 0%, / 0”, for the k, 
condition is given by: 


E (6) 


The in situ o stresses were determined from 
a geostatic stress field. Accordingly, 0. was 
computed using eq. (6) and the Poisson ratios, 
v, given in table 1. The values of c' and q” (ob- 


tained from equations (4) and (5), respectively), 


y 
3max 


as well as other strength parameters, computed 
for each material in the range 0<0%<0G”,... are 


given in table 2. 
Rock blocks contact properties 


The numerical method used for the analysis 
has a soft contact formulation, which uses 
normal, E and shear, k S stiffnesses and a 
yield criterion to simulate joint behavior. The 
Mohr-Coulomb joint model was used given 
its simplicity and its proven applicability 
to complex friction problems. The friction 


resistance T' (as a function of o”) between rocks 
blocks considered was T' = 0.85 o”, for o”, < 200 
MPa and Y = 0.50 + 0.60 o”, for o”, > 200 MPa 
(Byerlee, 1978). It is worth mentioning that the 
maximum principal stress value computed in 
the numerical model was much lower than 
200 MPa which in accordance with Byerlee 
(1978) would yield a friction coefficient of 
0.85 (p” approximately equal to 40%). Given the 
many uncertainties involved in the rock-rock 
friction problem, it was judged convenient to 
assume a conservative maximum value of 360. 
The minimum value was assumed 25% which 
corresponds to the minimum value obtained 
through the Hoek-Brown criterion (see table 
2). Table 3 shows the actual values used for 
each material contact. It is important to stress 
the fact that when two different materials 
were in contact, the lowest friction angle was 
considered. 

The magnitude of k, is computed through 
k, assuming an elastic relation between them 
(Kulhawy, 1975): 


k 
y k, =2k,(1+v) (7) 


Where v the joint material Poisson's ratio, 
and k, is the shear stiffnesses included in table 


(3). 


Table 2. Strength parameters of rock material. 


Material c' (MPa) 04 O, (MPa) 
1. Dacite V, = 3.8 km/s 2.84 49 8.60 
2. Dacite V, =3.2 km/s 4.36 49 14.13 
3. Dacite V, =2.7 km/s 2.41 47 5.42 
4. Fractured breccia VA =2.5km/s 1.68 44 3.33 
5. Uncompressed dacite 1.01 41 1.13 
6. Uncompressed breccia 0.91 39 0.64 
7. Dacite-rhyolite v, =1.7km/s 0.34 34 0.39 
8. Dyke (andesite) 0.03 25 0,14 
9. Vertedor fault 0.15 26 0.22 
10. Vitreous tuff 0.28 33 0.28 
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Table 3. Contact properties. 


Material Ñ eS la AS ES 
1. Dacite V, = 3.8 km/s 36 1:79 0.004 1168 433 
2. Dacite V, = 3.2 km/s 36 1.62 0.004 1093 405 
3. Dacite V, =2.7 km/s 36 1.80 0.005 972 360 
4. Fractured breccia V, =2.5 km /s 30 1.94 0.006 872 329 
5. Uncompressed dacite 30 2.07 0.006 934 346 
6. Uncompressed breccia 30 241) 0.007 820 304 
7. Dacite-rhyolite V, = 1.7 km/s 30 2.22 0.007 854 316 
8. Dyke (andesite) 25; pESS) 0.015 420 155 
9. Vertedor fault 25 2.31 0.015 416 154 
10 Vitreous tuff 27 2.28 0.015 410 [52 


Numerical model 


The numerical model was properly designed so 
that its boundary conditions would not affect 
the results of the excavation process (figure 
3). The model consisted of 225 955 tetrahedral 
solid elements of constant strain, and 88 655 
nodes. The model was built from 11 geological 
sections roughly spaced every 50 m along the 
spillway channel (see figure 1). From elevation 
544 to 624 amsl (above mean sea level), the left 
slope of the excavation was modeled with an 
inclination of 0.88:1 (H:V); from 624 to 700 ams] 


the slope inclination was 1.40:1; and from 700 
ams] to the top of the slope, the inclination was 
1.33:1. The right slope of the spillway channel 
had a constant inclination along its height 
(elevation 544 to 580 amsl) of 0.54:1. 

The numerical model is based on a time 
marching solving procedure finite difference 
which allows to explicitly simulating the 
contacts and discontinuities between blocks. 
The interaction between blocks is accounted for 
through stress-displacement relationships for 
each contact. Internal stresses in each block are 
computed through the finite difference method. 


Excavated material 


225 955 elements 
88 655 nodes 


Elev. 380 


(a) 


17 m wide berm 
Elev. 700 


(b) 


Figure 3. 3D geometry of the numerical model. 
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The most important geological discontinuities 
and faults were included in the model, while 
the minor rock fractures interacting along the 
major discontinuities were taken into account 
considering the fractured rocks as a rock mass 
(Hoek é€z Brown, 1997). 


Site initial conditions 


The results of the static and seismic (OBE and 
MCE seismic conditions) calculations of the 
site initial conditions showed no indication of 
any instability, which agrees with the historical 
stability of the rock mass. 


Excavation process simulation 


The excavation process from elevation 827 
to 544 amsl was simulated in four stages. 
Displacements and factors of safety were 
computed for each stage to evaluate any 
potential slope sliding during the excavation 
process. The results indicate that after the 
excavation was completed, the displacement 
field was mainly due to expansion of the 
rocks caused by stress relaxation, as shown 
in figure 4(a) through 4(d). The analysis show 


that the excavation process most likely will 
not induce massive rock falls. However, in 
figure 5(a) through 5(d) it is seen that the safety 
factor in some areas of the slope indicates 
that instabilities are likely to occur during 
the construction process. Accordingly, it was 
considered that important amounts of debris 
could fall down hill, jeopardizing the stable 
construction of the spillway slopes. 


Effect of excavation process in the rock mass 
and its implication in the factor of safety 


Due to the excavation process the rock 
mass suffers some damage causing strength 
reductions mainly near to the slope surface. 
Accordingly, a damage factor of 0.7 was 
considered, which corresponds to a good 
quality excavation process (Brady é Brown, 
2004). Figure 5d shows the factor of safety 
distribution in the maximum slope section for 
this damaged condition. It is noticeable that a 
zone of low factor of safety develops at the top 
of the slope. These results suggest that potential 
rock instabilities are likely to unfold during the 
construction process, thus indicating that the 
excavation should be properly monitored to 
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Displacement vectors 


(magnitude in meters) 
min max 
2.000E-02 2.400E-02 
1.600E-02 2.000E-02 
1.200E-02 1.600E-02 
8.000E-03 1.200E-02 


4.000E-03 8.000E-03 


0.000E+00 4.000E-03 


a) End of excavation stage 1 
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1.500E-02 2.000E-02 
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b) End of excavation stage 2 


Displacement vectors 
(magnitude in meters) 
min max 


2.500E-02 3.000E-02 
2.000E-02 2.500E-02 
1.500E-02 2.000E-02 
1.000E-02 1.500E-02 
5.000E-03 1.000E-02 
0.000E+00 5.000E-03 


c) End of excavation stage 3 


d) End of excavation stage 4 


Figure 4. Resultant displacement field due to excavation process. 
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Factor of safety 
max min 
1.000E+00 5.000E-01 
1.500E+00 1.000E+00 
2.000E+00 1.500E+00 
2.000E+00 2.000E+00 
2.000E+00 2.000E+00 
.. 2.000E+00 2.000E+00 


a) End of excavation stage 1 


Factor of safety 
max min 
1.000E+00 5.000E-01 
1.500E+00 1.000E+00 
2.000E+00 1.500E+00 
2.000E+00 2.000E+00 
2.000E+00 2.000E+00 
2.000E+00 2.000E+00 


c) End of excavation stage 3 


Factor of safety 
max min 
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2.000E+00 2.000E+00 
2.000E+00 2.000E+00 
2.000E+00 2.000E+00 


b) End of excavation stage 2 


Factor of safety 
max min 


1.000E+00 5.000E-01 
1.500E+00 1.000E+00 
2.000E+00 1.500E+00 
2.000E+00 2.000E+00 
2.000E+00 2.000E+00 
2.000E+00 2.000E+00 


d) End of excavation stage 4 


Figure 5. Evolution of safety factor due to excavation process. 


detect any sign of rock movements and proceed 
to stabilize the area as required. 


Seismic analysis 


Despite the effect of transient loads can 
affect noticeably the behavior of rock masses, 
the physics of the phenomenon is poorly 
understood (Jing, 2003). Henceforth, linear 
elastic rock models are widely accepted for 
analyzing seismic rock mass responses. The 
dynamic response of the slope was computed 


for both the OBE and the MCE according to the 
design criteria used for La Yesca dam (ICOLD, 
2010). Synthetic earthquakes were developed to 
obtain the OBE and MCE time series. The input 
motions computed were representative of the 
Mexican subduction (Cocos plate) earthquake- 
generating mechanism. 

The seismic response of the slope at the 
end of the excavation for the OBE earthquake 
is shown in figure 6, where it is seen that 
displacements up to 40 cm of were developed 
at the top of the slope. This indicates that 


Displacement vectors 
(magnitude in meters) 


3.750E-01 4.500E-01 
3.000E-01 3.750E-01 
6.250E-02 7.500E-02 
5.000E-02 6.250E-02 
3.750E-02 5.000E-02 
2.500E-02 3.750E-02 
1.250E-02 2.500E-02 
1.200E-02 1.600E-02 
8.000E-03 1.200E-02 
4.000E-03 8.000E-03 
0.000E+00 4.000E-03 


Figure 6. Resultant displacements at the end of the excavation. 
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instabilities can occur over that area, which 
in case of rock debris sliding towards the 
spillway channel, could lead to its obstruction 
or to a malfunction of it, endangering the dam 
itself. Also, the results show that a potentially 
unstable block can develop at the toe of the 
slope. 

The MCE is a highly severe seismic event. 
The resultant slope displacements by action 
of this earthquake increase dramatically 
with respect to the OBE earthquake. At the 
top of the slope the displacements indicate 
the failure in this zone of the rock cut (figure 
7). Displacements at the toe of the slope also 
increase strongly. Displacement magnitudes at 
the slope toe reach about 40 cm, thus indicating 
possible block detachments under the action of 
the maximum credible earthquake. That could 
lead to a major slope slide towards the spillway 
channel. 


Corrective actions 


The analysis results show the potential 
instabilities of the spillway slopes and hence 
impairing its correct operation. Accordingly 
it was decided to modify both its initial lay 
out and geometrical design. To this end The 
Federal Commission of Electricity (owner of 
the project) implemented important changes 
observed when comparing figure 8 (original 
design) with figures 9 and 10 (final design). The 


main modification on the spillway lay out was 
changing the orientation (with respect to the 
original) in 28? to the southeast which reduced 
the slope maximum height by 15 m. This had a 
benefic effect in decreasing excavation volumes 
and lessening the hazard of falling rock blocks 
into the spillway channel that could likely 
impair its proper operation or even clog it, 
putting at a high risk the dam safety. Regarding 
the modification of the original spillway 
consisted in dividing the one channel section 
into a three channel section, which in addition 
to increasing the overall stability of the spillway 
slopes, the material volume excavation was 
decreased. Picturel0 depict a view of the final 
configuration of the La Yesca dam spillway. 


Conclusions 
+ From the analysis performed, it can be 


concluded that modeling 
is an important tool in the assessment 


numerical 


of geotechnical projects. Results from 
such analysis allow forestalling possible 
problems that can arise during project 
execution. Also, numerical modeling is a 
helpful tool for selecting the zones where 
detailed monitoring should be carried out. 
e However, the numerical models, as itis well 
known, should be fed with realistic material 
properties obtained from laboratory and / 
or field tests in order to obtain meaningful 


Displacement vectors 
(magnitude in meters) 


3.000E+00 3.750E+00 
2.250E+00 3.000E+00 
5.000E-01 6.000E-01 
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3.000E-01 4.000E-01 
2.000E-01 3.000E-01 
1.000E-01 2.000E-01 
0.000E+00 1.000E-01 


Figure 7. Resultant slope displacements by action of MCE. 
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Figure 9. Final spillway design. 


Figure 10. Final spillway conformation (Modified www.cfe.gob.mx). 
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results. When material properties are not 
known properly, parametric studies should 
be carried outto evaluate theinfluence ofthe 
various parameters. Itis worth mentioning 
that the spillway, along its new trace and 
the geometry changes mentioned, was 
executed without any major problems. It is 
worth mentioning that numerical analysis 
allow determination of “hidden” problems, 
which with other procedures may not be 
detected. Accordingly, numerical methods 
are powerful tools that provide clues for 
taking educated design decisions and 
consequently improve structural security 
and help lowering risk levels. 

+ The implemented changes were very 
important, and the visualization of the 
potential problems only was possible with 
the 3D numerical study. 
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